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Abstract 


Maximum likelihood (ML) and Bayesian inference (BI) analyses using two mitochondrial (16S and cyt 5) and two nuclear (CMOS 
and RAG1) genes and 103 specimens recovered the first phylogenies of all 23 extant species of Goniurosaurus. The analyses strongly 
supported the recognition of four monophyletic species groups with identical inter-specific relationships within the kuroiwae, lichten- 
felderi, and yingdeensis groups but discordant topologies at some nodes within the /uii group. Both analyses recovered a polyphyletic 
G. luii with respect to G. kadoorieorum, and owing to the lack of diagnostic characters in the latter, it is considered a junior synonym 
of G. luii. A stochastic character mapping analysis of karst versus non-karst habitat preference suggested that karstic landscapes may 
have played a major role in the evolution and diversification of Goniurosaurus. A karst habitat preference is marginally supported as 
the most probable ancestral condition for Goniurosaurus as well as for the kuroiwae, luii, and yingdeensis groups. However, a non- 
karst habitat preference is marginally supported as the most probable ancestral condition for the /ichtenfelderi group. Multivariate 
and univariate ecomorphological analyses of the karst-adapted G. catbaensis, G. huuliensis, and G. luii of the luii group and the gran- 
ite-stream-adapted G. lichtenfelderi of the lichtenfelderi group demonstrated that their markedly statistically different body shapes 
may be an adaptive response that contributes to habitat partitioning in areas of northern Vietnam where they are nearly sympatric. 
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Introduction 


Eublepharid geckos of the genus Goniurosaurus Barb- 
our, 1908 comprise 23 saxicolous specialists (Uetz et al. 
2021) that extend from the Ryukyu Archipelago in Japan, 
southward through East Asia to northern Vietnam. Goni- 
urosaurus is a well-defined monophyletic group (Grismer 
1988) comprised of four monophyletic species groups: 
the kuroiwae group containing six species endemic to 
the Ryukyu Archipelago, Japan; the lichtenfelderi group 
with five species from insular and mainland China and 
northern Vietnam; the /uii group with eight species from 
northern Vietnam, some of its offshore islands, and south- 
ern China; and the yingdeensis group consisting of four 
species from southern China (Kurita et al. 2008; Nguyen 
et al. 2009; Nguyen 2011; Wang et al. 2013; Honda and 
Ota 2017; Liang et al. 2018; Qi et al. 2020a, 2020b; Zhu 
et al. 2020a, 2020b). Apart from these species, Goniuro- 
saurus sinensis Zhou, Peng, Huo and Yuan, 2019 is likely 
a junior synonym of another species from Hainan Island, 
China and not included herein (Qi et al. in progress). Phy- 
logenetic relationships within Goniurosaurus have never 
been strongly supported nor consistent among different 
studies (e.g. Wang et al. 2013; Liang et al. 2018; Qi et al. 
2020a, 2020b; Zhu et al. 2020a, 2020b). This protract- 
ed state of discordance results, in part, from researchers 
focusing on different species groups as opposed to the 
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entire genus, as well as using different genes or different 
combinations of genes with varying combinations of in- 
group and outgroup species—all variables that bear sig- 
nificantly on tree construction (Wiens 1998; Zwickl et al. 
2002; Heath et al. 2008; Wiens and Morrill 2011; Wain- 
wright and Price 2016). The most commonly used genet- 
ic markers have been the mitochondrial genes 12S and 
16S rRNA and cytochrome b (cyt 5). Liang et al. (2018) 
were the first to address the challenges of properly align- 
ing rRNA (Pyron et al. 2013) and constructed a well-sup- 
ported mito-nuclear data set using 16S, cyt b, and the nu- 
clear genes oocyte maturation factor MOS (CMOS), and 
recombination activating | (RAG1). Zhu et al. (2020b) 
also used this mito-nuclear combination, but examined 
only relationships within the lichtenfelderi group. 

In an effort to continue building a more global under- 
standing of the phylogenetic relationships within Goni- 
urosaurus, We expanded the mito-nuclear data set of Li- 
ang et al. (2018) to include 103 individuals as opposed 
to 31 and 23 as opposed to 17 species, which for the first 
time, includes all extant species of the genus (Table 1). We 
used this phylogeny in a stochastic character state map- 
ping (SCM) analysis (Revell 2012) of habitat preference 
to explore the role karstic landscapes may have played in 
the evolution and diversification of Goniurosaurus and if 
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Figure 1. Distribution of the four recorded species of Goniurosaurus in Vietnam and China used in the ecomorphological analysis. 
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Granite-stream-adapted species (G. lichtenfelderi), B1 — B3. Karst habitats of the /uii group, B4. Granite-stream habitat of the 


lichtenfelderi group. Photos by Hai Ngoc Ngo. 
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Table 1. Species and GenBank accession numbers of the sequenced specimens used herein. 


Species/Specimen 16s cytb CMOS RAGI 
Goniurosaurus araneus AB308460 

G. araneus ECNU-V0008 MT533259 

G. araneus JFBM15830 HQ426537 HQ426286 
G. bawanglingensis BL-RBZ-021 MH247190 MH247201 MH247212 MH247223 
G. bawanglingensis BL-RBZ-022 MH247191 MH247202 MH247213 MH247224 
G. bawanglingensis BL-RBZ-023 MH247192 MH247203 MH247214 MH247225 
G. bawanglingensis BL-RBZ-024 MH247193 MH247204 MH247215 MH247226 
G. bawanglingensis SYS 1002162 MT995758 MT995773 

G. catbaensis G33 MW741550 MwW650944 

G. catbaensis G34 MW741551 MwW650945 

G. catbaensis G35 MW650946 

G. catbaensis MHNG 2699.49 EU499389 

G. gezhi ECNU-V0038 MT533260 

G. gezhi ECNU-V0040 MT533261 

G. gezhi ECNU-V0042 MT533262 

G. gezhi ECNU-V0046 MT533263 

G. gezhi ECNU-V0047 MT533264 

G. gollum SYS 1002420 MT995784 MT995787 MW727559 MW727594 
G. gollum SYS 1002421 MT995785 MT995788 MW727560 MW727595 
G. gollum SYS 1002422 MT995786 MT995789 MW727561 MW727596 
G. hainanensis BL-RBZ-041 MH247194 MH247205 MH247216 MH247227 
G. hainanensis BL-RBZ-042 MH247195 MH247206 MH247217 MH247228 
G. hainanensis SYS 1000349 KC765080 

G. hainanensis JK1 AB308458 

G. huuliensis Gohu AB853453 AB853479 

G. huuliensis G21 MW650936 

G. huuliensis G23 MW650937 

G. huuliensis G24 MW650938 

G. kadoorieorum ECNU-V0058 MT533258 

G. kadoorieorum ECNU-V0060 MT533265 

G. kadoorieorum ECNU-V0061 MT533266 

G. kuroiwae Gokul Northern Okinawa AB853448 AB853473 

G. kuroiwae Goku2 Southern Okinawa AB853445 

G. kuroiwae Goor| Southern Okinawa AB853446 AB853467 

G. kwanghua ECNU-V0003 MK782788 MK782782 MK782776 MK782770 
G. kwanghua ECNU-V0004 MK782789 MK782783 MK782777 MK782771 
G. kwanghua ECNU-V0005 MK782790 MK782784 MK782778 MK782772 
G. kwangsiensis ECNU-V0009 MK782786 MK782780 MK782774 MK782768 
G. liboensis SYS 1000217 KC900230 

G. lichtenfelderi ECNU-V0007 MK782785 MK782779 MK782773 MK782767 
G. lichtenfelderi IEBR 3692 JF799756 

G. luii ECNU-V0012 MK782787 MK782781 MK782775 MK782769 
G. luii Golu2 EF081254 

G. luii Golu3 AB853452 AB853478 

G. luii SYSr 000255 KC765083 

G. luii SYSr 000256 KC765084 

G. luii ZFMK 87057 EU499391 

G. luii TG00795 HQ426287 
G. orientalis Goku3 AB853446 

G. orientalis Goor2 AB853443 AB853461 

G. orientalis Goor3 AB853462 

G. sengokui Gose1 AB853444 AB853463 

G. sengokui Gose2 AB853464 

G. sengokui KUZ 62087 HQ876393 

G. splendens Gosp1 AB853451 AB853477 

G. splendens Gosp2 AB853449 

G. splendens Gosp3 AB853450 
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Species/Specimen 16s cytb CMOS RAGI 

G. toyamai Gotol AB853447 AB853468 

G. toyamai Goto2 AB853469 

G. varius SYS 1002331 MT995754 MT995769 MW727556 MW727590 
G. varius SYS 1002332 MT995755 MT995770 

G. varius SYS 1002333 MT995753 MT995768 

G. varius SYS 1002362 MT995756 MT995771 MW727557 MW727592 
G. varius SYS 1002363 MT995757 MT995772 MW727558 MW727593 
G. varius SYS 1002485 MW721828 MW727532 MW727562 MW727597 
G. varius SYS 1002486 MW721829 MW727533 MW727563 MW727598 
G. yamashinae Goyal AB853442 AB853460 

G. yamashinae Goya2 AB853441 AB853459 

G. yamashinae Goya3 AB853458 

G. yingdeensis Field number DYAO1 MW721830 MW727534 MW727574 MW727605 
G. yingdeensis Field number DYA02 MW721831 MW727535 MW727575 MW727606 
G. yingdeensis Field number HSO1 MW721832 MW727536 MW727576 MW727607 
G. yingdeensis Field number HS02 MW721833 MW727537 MW727577 MW727608 
G. yingdeensis Field number LTO1 MW721834 MW727538 MW727580 MW727611 
G. yingdeensis Field number LT02 MW721835 MW727539 MW727581 MW727612 
G. yingdeensis SYS 1000549 KC765082 

G. yingdeensis SYS 1000550 KC900231 

G. yingdeensis SYS 1001271 MT995759 MT995774 MW727547 

G. yingdeensis SYS 1001272 MT995760 MT995775 MW727548 

G. yingdeensis SYS 1001493 MT995761 MT995776 MW727551 

G. yingdeensis SYS 10002115 MT995762 MT995777 

G. zhelongi Field number HW01 MW721838 MW727540 MW727578 MW727609 
G. zhelongi Field number HW02 MW721839 MW727541 MW727579 MW727610 
G. zhelongi Field number MDAO1 MW721836 MW727542 MW727582 MW727613 
G. zhelongi Field number MDA02 MW721837 MW727543 MW727583 MW727614 
G. zhelongi Field number TZO1 MW721840 MW727544 MW727584 MW727615 
G. zhelongi Field number TZ02 MW721841 MW727545 MW727585 MW727616 
G. zhelongi SYS 1000816 KJ423105 MT995778 MW727570 

G. zhelongi SYS 1001491 MT995763 MT995779 MW727549 

G. zhelongi SYS 1001492 MT995764 MT995780 MW727550 

G. zhelongi SYS 1002108 MT995765 MT995781 

G. zhoui BL-RBZ-001 MH247196 MH247207 MH247218 MH247229 
G. zhoui BL-RBZ-004 MH247197 MH247208 MH247219 MH247230 
G. zhoui BL-RBZ-006 MH247198 MH247209 MH247220 MH247231 
G. zhoui BL-RBZ-007 MH247199 MH247210 MH247221 MH247232 
G. zhoui BL-RBZ-008 MH247200 

G. zhoui SYS 1002213 MT995766 MT995782 MW727553 

G. zhoui SYS 1002214 MT995767 MT995783 MW727554 

Eublepharis macularius AB853454 AB853480 


habitat preference coevolved with ecomorphology in near 
sympatric species of the /uii and lichtenfelderi groups in 
Vietnam (Ngo et al. 2021; Figs. 1, 2). 


Materials and methods 


Genetic data and phylogenetic 
analyses 


Genomic DNA was extracted from muscle tissue sam- 
ples, using a DNA extraction kit from Tiangen Biotech 
(Beijing) Co., Ltd. Primers used for 16S were rl16S-5L 


(5’°- GGIMMYGCCTGCCCAGTG -3’) and 16Sbr-H 
(5’- CCGGTCTGAACTCAGATCACGT-3’) (Palumbi 
et al. 1991), for cyt b the primers were L14731 (5’- TG- 
GTCTGAAAAACCATTGTTG-3’) (Honda et al. 2014) 
and H15149m (5’- GCMCCTCAGAAKGATATTTGY- 
CCTCA-3’) (Chambers and MacAvoy 1999), for CMOS 
the primers were FU-F (5’- TTTGGTTCKGTCTACAA- 
GGCTAC -3’) and FU-R (5’- AGGGAACATCCAAAG- 
TCTCCAAT -3’) (Gamble et al., 2008) , and for RAGI 
the primers were R13 (5’- TCTGAATGGAAATTCAAG- 
CTGTT -3’) and R18 (5’- GATGCTGCCTCGGTCGG- 
CCACCTTT -3’) (Groth and Barrowclough 1999). The 
PCR procedure was performed with an initial denatur- 
ation at 94 °C for 5 min, 35 cycles of 94 °C for 30 s, 
55 °C for 30 s and 72 °C for | min, followed by a final 
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extension at 72 °C for 10 min (Liang et al. 2018). PCR 
products were purified with spin columns and then se- 
quenced with forward primers using BigDye Terminator 
Cycle Sequencing Kit as per the guidelines on an ABI 
Prism 3730 automated DNA sequencer by Shanghai Ma- 
jorbio Bio-pharm Technology Co., Ltd. 

We constructed Maximum Likelihood (ML), Bayes- 
ian Inference (BI), and Bayesian Evolutionary Analysis 
by Sampling Trees (BEAST) phylogenetic trees using a 
concatenated data set composed of 3070 base pairs (bp) 
of the mitochondrial genes, 16S (633 bp) and cyt b (1075 
bp), and the nuclear genes, CMOS (472 bp) and RAGI 
(890), from 103 specimens of 23 species of Goniurosau- 
rus with varying degrees of sequence coverage across the 
samples (Table 1). Concatenation followed the compari- 
son of separate gene trees to confirm there were no ma- 
jor discordances. One species, Eublepharis macularius, 
served as an outgroup (Grismer 1988; Jonniaux and Ku- 
mazawa 2008) to root the trees. Sequence data and Gen- 
Bank accession numbers are listed in Table 1. 

A Maximum likelihood (ML) analysis partitioned 
by gene was implemented using the IQ-TREE web- 
server (Nguyen et al. 2015; Trifinopoulos et al. 2016) 
preceded by the selection of substitution models using 
TIM2+F+I+G4 for 16S and cyt b and HK Y+F for CMOS 
and RAGI. To avoid over parameterization, protein cod- 
ing genes were not partitioned by codon. One-thousand 
bootstrap pseudoreplicates via the ultrafast bootstrap 
(UFB: Hoang et al. 2018) approximation algorithm were 
employed, and nodes having UFB values of 95 and above 
were considered strongly supported (Minh et al. 2013). 
We considered nodes with values of 90-94 as well-sup- 
ported. A Bayesian inference (BI) analysis was carried 
out in MrBayes 3.2.3. (Ronquist et al. 2012) on XSEDE 
using the CIPRES Science Gateway (Cyberinfrastructure 
for Phylogenetic Research; Miller et al. 2010) employing 
default priors and models of evolution that most close- 
ly approximated those selected by the BIC and used in 
the ML analysis. Two independent Markov chain Mon- 
te Carlo (MCMC) analyses for each data set were per- 
formed—each with four chains, three hot and one cold. 
The MCMC simulations ran for 100 million generations. 
Trees were sampled every 10,000 generations, and the 
first 10% of the trees from each run from each data set 
were discarded as burn-in. The parameter files from 
both runs were checked in Tracer v1.6 (Rambaut et al. 
2014) to ensure that convergence and stationarity of all 
parameters had effective sample sizes (ESS) well-above 
above 200. Post-burn-in sampled trees from each respec- 
tive run were combined and 50% majority-rule consensus 
trees were constructed. Nodes with Bayesian posterior 
probabilities (BPP) of 0.95 and above were considered 
highly supported (Huelsenbeck et al. 2001; Wilcox et al. 
2002). We considered nodes with values of 0.90—0.94 as 
well-supported. 

An input file was constructed in Bayesian Evolution- 
ary Analysis Utility (BEAUti) v. 2.4.6 using a relaxed 
lognormal clock with unlinked site models, linked trees 
and clock models, and a Yule prior and run in BEAST on 
CIPRES (Cyberinfrastructure for Phylogenetic Research; 


Miller et al. 2010). bModelTest was used to numerically 
integrate over the uncertainty of substitution models of 
each gene while simultaneously estimating phylogeny us- 
ing Markov chain Monte Carlo (MCMC). MCMC chains 
were run for 100,000,000 generations and logged every 
10,000 generations. The BEAST log file was visualized 
in Tracer v. 1.6.0 (Rambaut et al. 2014) to ensure effec- 
tive sample sizes (ESS) were well-above 200 for all pa- 
rameters. A Maximum clade credibility tree using mean 
heights at the nodes was generated using TreeAnnotator 
v.1.8.0 (Rambaut and Drummond 2014) with a burn-in 
of 1000 trees (10%). Nodes with BPPs of 0.95 and above 
were considered strongly supported (Huelsenbeck et al. 
2001; Wilcox et al. 2002). We considered nodes with val- 
ues of 0.90—0.94 as well-supported. 


Ancestral state reconstruction 


The BEAST tree was converted to newick format and 
pruned using the drop.tip() command (Paradis and 
Schliep 2018) in the R package ape [v.3.4.3] to include 
only the earliest diverged individual of each species. 
Habitat preference (karst or non-karst; see below) was 
mapped onto the tree using SCM implemented in the R 
package Phytools (Revell 2012) in order to derive prob- 
ability estimates of the ancestral states at each node. A 
transition rate matrix was identified that best fit the data 
by comparing the corrected Akaike Information Criterion 
(AICc) values in the R package ape (Paradis and Schilep 
2018). Three transition rate models were considered: a 
2-parameter model having different rates for every tran- 
sition type (the ARD model); a single-parameter model 
with equal forward and reverse rates between states (the 
symmetrical rates SYM model); and a single rate param- 
eter model that assumes equal rates among all transitions 
(ER). Lastly, an MCMC approach was used to sample the 
most probable 1000 trait histories from the posterior us- 
ing the make.simmap() command and then summarized 
them using the suwmmary() command. 

The coding of habitat preference for each species was 
determined from the literature and field observations of 
the authors (Table 2). A species’ habitat preference was 
coded as “karst” if it had a strong association with karstic 
habitats. Many such species may range into forested ar- 
eas or even other rock types (granite or volcanic) if near- 
by but their densities are much greater in karstic areas. 
A species was coded as “non-karst” if detected only in 
forested areas or forested areas with other rock types (e.g. 
granite). These species never show any strong preference 
for karstic microhabitats even if such habitats exist within 
their range. 


Morphological data and analyses 


An ecomorphological analysis was conducted using 
four of the five recorded species from Vietnam (Grismer 
et al. 1999; Vu et al. 2006; Orlov et al. 2008; Ziegler et 
al. 2008; Nguyen et al. 2009; Nguyen 2011; Orlov et al. 
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Table 2. Habitat preferences of the species of Goniurosaurus. 


Species 1° habitat 2° habitat Source 

kuroiwae group 

G. splendens karst forest Nakamura and Ueno (1963), H. Ota pers. comm., L. Grismer pers. obs. 
G. toyamai forest H. Ota pers. comm. 

G. kuroiwae North forest H. Ota pers. comm., L. Grismer pers. obs., 
G. kuroiwae South karst forest Nakamura and Ueno (1963), H. Ota pers. comm., L. Grismer pers. obs. 
G. yamashinae karst forest H. Ota pers. comm., L. Grismer pers. obs., 
G. sengokui karst forest Werner et al. (2004), H. Ota pers. comm. 
G. orientalis karst H. Ota pers. comm. 

yingdeensis group 

G. gollum karst Qi et al. (2020a) 

G. yingdeensis karst granite Wang et al. (2010), S. Qi pers. obs. 

G. zhelongi karst granite S. Qi, pers. obs., Wang et al. (2014) 

G. varius karst Qi et al. (2020b) 

lichtenfelderi group 

G. bawanglingensis granite karst Grismer et al. (2002), Orlov et al. (2008) 
G. zhoui karst granite Zhou et al. (2018), S. Qi pers. obs. 

G. kwanghua karst granite Zhu et al. (2020) 

G. lichtenfelderi granite Orlov et al. (2008) 

G. hainanensis granite volcanic S. Qi pers. obs., L. Grismer pers. obs. 

luii group 

G. catbaensis karst Ziegler et al. (2008), Ngo et al. (2019a) 

G. gezhi karst Zhu et al. (2020) 

G. araneus karst Grismer et al. (1999) 

G. kadoorieorum karst Yang and Chan (2015) 

G. huuliensis karst Orlov et al. (2008) 

G. luii karst Grismer et al. (1999), Vu et al. (2006) 

G. liboensis karst Wang et al. (2013) 

G. kwangsiensis karst Yang and Chan (2015) 


2020) for which there existed a substantially large mor- 
phometric data set (Ngo et al. 2021): the karst-adapted G. 
catbaensis Ziegler, Nguyen, Schmitz, Stenke, and Résler, 
2008, G. huuliensis Orlov, Ryabov, Nguyen, Nguyen, and 
Ho, 2008, and G. /uii Grismer, Viets, and Boyle, 1999 
of the /uii group and the granite stream-adapted G. Jicht- 
enfelderi (Mocquard, 1897) of the lichtenfelderi group 
(Figs. 1, 2). A total of 486 live individuals and 54 museum 
specimens of four species were examined for morpholog- 
ical data, comprising 194 individuals of G. catbaensis (21 
juveniles, 93 females, and 80 males), 80 individuals of G. 
huuliensis (two juveniles, 46 females, and 32 males), and 
88 individuals of G. /uii (11 juveniles, 43 females, and 34 
males) of the /uii species group and 178 individuals of G. 
lichtenfelderi (14 juveniles, 72 females, and 92 males) of 
the lichtenfelderi group. 

Measurements were taken with dial calipers to the 
nearest 0.1 mm on the right side of each individual. Ab- 
breviations are as follows: snout-vent length (SVL), from 
tip of snout to vent; axilla to groin length (AG), from 
posterior edge of forelimb insertion to anterior edge of 
hind limb insertion; maximum body width (BW), greatest 
width of torso, taken at level of midbody; maximum body 
height (BH), from dorsal surface of body to belly; inter- 
narial distance (ID), distance between nares; head length 
(HL), from the tip of snout to posterior edge of occiput; 
maximum head width (HW); cheek height (CH), from 


posterior edge of labial to top of head at parietal region; 
interorbital distance (IO), distance between posteriormost 
points of eyes; diameter of auditory meatus (AD); snout 
to eye distance (SL), measured from tip of snout to an- 
teriormost point of eye; diameter of eye (ED), greatest 
diameter of eye; eye to ear distance (EE), from posterior 
margin of eye to posterior margin of ear; forelimb length 
(FLL), from axilla to the tip of the fourth finger; hind 
limb length (HLL), from groin to the tip of the fourth toe. 
To remove potential effects of allometry, size was adjust- 
ed using the following equation: X,,=log(X)-B[log(SVL 
)-log(SVL jean) ], Where X,4=adjusted value; X=measured 
value; B=unstandardized regression coefficient for each 
population; and SVL,,,.., overall average SVL of all pop- 
ulations (Thorpe 1975, 1983; Turan 1999; Lleonart et al. 
2000)—accessable in the R package GroupStruct (avail- 
able at https://github.com/chankinonn/GroupStruct). The 
morphometrics of each species were adjusted separately 
and then concatenated so as not to conflate intra- with 
interspecific variation (Reists 1986). All data were then 
scaled to their standard deviation to insure they were an- 
alyzed on the basis of correlation and not covariance and 
were log-transformed to insure they were normally dis- 
tributed. 

An analysis of variance (ANOVA) was performed on 
a data set coded for species to search for the presence 
of statistically significant mean differences (p < 0.05) 
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Figure 3. Mito-nuclear maximum likelihood topology with ultrafast bootstrap values (UFB) and Bayesian posterior probabilities 
(BPP) at the nodes. All species except Goniurosaurus luii had strong nodal support (100/1.00) for their monophyly. The inset in the 
uii species group is a section of the BI analysis showing the non-monophyly of G. /uii with respect to G. kadoorieorum. Colored 
species are those used in the ecomorphological analyses. 
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Figure 4. Stochastic character map of karst and non-karst habitat preferences on a BEAST topology. Percentages at the nodes are 


the probability of the ancestral habitat preference. Numbers at the nodes are BPP values. Colored species are those used in the eco- 


morphological analyses. 


among characters across the selected subset of species in 
the /uii and lichtenfelderi groups. Character means bear- 
ing statistical differences among species were subjected 
to a TukeyHSD test to ascertain which species pairs dif- 
fered significantly from each other for those particular 
characters. A Student ¢-test was also performed on a sec- 
ond data set coded for only habitat preference (karst ver- 
sus non-karst) to search for the presence of statistically 
significant mean differences (p < 0.05) among the same 
subsets of species coded for habitat. Violin plots with in- 
serted boxplots were generated in order to visualize the 
range, frequency, mean, 50% quartile, and degree of dif- 
ferences between the dependent variables for both data 
sets bearing statistically different mean values. 

The morphospatial clustering of the two separate data 
sets (species and habitat preference) were visualized us- 
ing principal component analysis (PCA) along the ordi- 
nation of the first two principal components (PC) using 
the Adegenet package in R (Jombart et al. 2010) and 
implemented by the prcomp() command. The data were 
log-transformed prior to analysis in order to normalize 
their distribution so as to ensure characters with very 
large or very low values could not over-leverage the re- 
sults owing to intervariable nonlinearity. All statistical 
analyses were performed using R.3.1.2 (R Core Team 
2018). 


Results 


Phylogenetic relationships 


The ML, BI, and BEAST analyses recovered strong nod- 
al support (UFB 98—100/BPP 1.00) for the monophyly 
of all four species groups with the kuroiwae group being 
the strongly supported (100/1.00) sister group to the re- 
maining three groups (Fig. 3). The ML analysis weakly 
recovered (88) the lichtenfelderi and yingdeensis groups 
as sister lineages but the BEAST analysis recovered this 
relationship with strong support (1.00; Fig. 4). The BI 
analysis recovered the /uii and yingdeensis groups as sis- 
ter lineages, although the support is so low (0.51), the 
three groups effectively form a polytomy. The ML and 
BI analyses recovered the identical inter-specific relation- 
ships within the species groups but discordant relation- 
ships with the BEAST analysis regarding the /uii group. 
The ML and BI analyses recovered a poorly supported 
(G. catbaensis (G. araneus, G. gezhi)) clade but the 
BEAST analysis recovered G. catbaensis as the strongly 
supported (0.99) sister species to the remainder of the /uii 
group species (Figs. 3, 4, respectively). All analyses re- 
covered a polyphyletic Goniurosaurus luii with respect to 
G. kadoorieorum (not shown in the pruned tree of Fig. 4). 
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Table 3. Difference, lower and upper ranges, and adjusted p values of statistically significant mean differences between species pairs 
for each character based on ANOVA and subsequent TukeyHSD analyses. 


axilla-groin (AG) difference lower range upper range p adjusted 
huuliensis-catbaensis 0.095633573 0.072292064 0.118975082 3.37E-10 
lichtenfelderi-catbaensis —0.075967576 —0.094200659 —0.057734493 3.37E-10 
luii-catbaensis —0.029676257 —0.05225407 —0.007098445 0.004206065 
lichtenfelderi-huuliensis —0.171601149 —0.195246969 —0.14795533 3.37E-10 
luii-huuliensis —0.12530983 —0.152447204 —0.098172456 3.37E-10 
luii-lichtenfelderi 0.046291319 0.023399042 0.069183596 1.60E-06 
body width (BW) 

huuliensis-catbaensis 0.106179184 0.066412755 0.145945613 4.37E-10 
lichtenfelderi-huuliensis —0.0833 14667 —0.123599543 —0.043029791 8.67E-07 
luii-huuliensis —0.121548632 —0.167781995 —0.0753 15269 5.34E-10 
body height (BH) 

huuliensis-catbaensis 0.094637915 0.046801122 0.142474708 2.84E-06 
lichtenfelderi-huuliensis —0.07317021 —0.121630666 —0.024709755 0.000647886 
luii-huuliensis —0.125013692 —0.180629845 —0.069397539 7.11E-08 
luii-lichtenfelderi —0.05 1843482 —0.098759605 —0.004927359 0.02359844 
internarial distance (ID) 

lichtenfelderi-catbaensis —0.082396274 —0.105742691 —0.059049857 3.37E-10 
luii-catbaensis —0.076638786 —0.105548379 —0.047729192 4.74E-10 
lichtenfelderi-huuliensis —0.096677642 —0.126954757 —0.066400527 3.37E-10 
luii-huuliensis —0.090920154 —0.125668004 —0.056172303 5.78E-10 
head length (HL) 

huuliensis-catbaensis 0.075818967 0.058801575 0.092836359 3.37E-10 
lichtenfelderi-catbaensis —0.162875997 —0.176169033 —0.149582961 3.37E-10 
luii-catbaensis —0.029098886 —0.045559496 —0.012638276 3.82E-05 
lichtenfelderi-huuliensis —0.238694964 —0.255934217 —0.221455712 3.37E-10 
luii-huuliensis —0.104917853 —0.124702663 —0.085133043 3.37E-10 
luii-lichtenfelderi 0.133777111 0.117087238 0.150466985 3.37E-10 
head width (HW) 

huuliensis-catbaensis 0.036775074 0.019138869 0.05441128 6.89E-07 
lichtenfelderi-catbaensis —0.158637831 —0.172414249 —0.144861413 3.37E-10 
luii-catbaensis —0.065002064 —0.082061241 —0.047942887 3.37E-10 
lichtenfelderi-huuliensis —0.195412905 —0.213279039 —0.177546772 3.37E-10 
luii-huuliensis —0.101777138 —0.122281395 —0.081272881 3.37E-10 
luii-lichtenfelderi 0.093635767 0.07633899 0.110932545 3.37E-10 
head height (HH) 

huuliensis-catbaensis 0.108413032 0.073172094 0.143653969 3.37E-10 
lichtenfelderi-catbaensis —0.032237965 —0.059766217 —0.004709713 0.014142658 
lichtenfelderi-huuliensis —0.140650997 —0.176351381 —0.104950613 3.37E-10 
luii-huuliensis —0.14002168 —0.180993602 —0.099049757 3.37E-10 
cheek height (CH) 

huuliensis-catbaensis 0.069595379 0.028246161 0.110944597 0.000101199 
lichtenfelderi-catbaensis —0.14807373 —0.180373429 —0.11577403 3.37E-10 
luii-catbaensis —0.057345812 —0.09734215 —0.017349473 0.001381599 
lichtenfelderi-huuliensis —0.217669109 —0.259557409 —0.175780808 3.37E-10 
luii-huuliensis —0.126941191 —0.175014741 —0.078867641 5.00E-10 
luii-lichtenfelderi 0.090727918 0.050174509 0.131281326 8.27E-08 
interorbital distance (ID) 

huuliensis-catbaensis 0.05692565 0.034698839 0.079152461 9.29E-10 
lichtenfelderi-catbaensis —0.211451872 —0.228814215 —0.194089528 3.37E-10 
luii-catbaensis —0.03636425 —0.057863835 —0.014864664 9.22E-05 
lichtenfelderi-huuliensis —0.268377522 —0.29089411 —0.245860933 3.37E-10 
luii-huuliensis —0.093289899 —0.1191313 —0.067448499 3.37E-10 
luii-lichtenfelderi 0.175087622 0.15328859 0.196886655 3.37E-10 
snout length (SL) 

huuliensis-catbaensis 0.108547374 0.083797596 0.133297152 3.37E-10 
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lichtenfelderi-catbaensis —0.20706889 —0.226402034 —0.187735746 3.37E-10 
lichtenfelderi-huuliensis —0.315616264 —0.340688712 —0.290543816 3.37E-10 
luii-huuliensis —0.125435452 —0.154210112 —0.096660792 3.37E-10 
luii-lichtenfelderi 0.190180812 0.16590737 0.214454254 3.37E-10 
ear diameter (ED) 

huuliensis-catbaensis —0.137477526 —0.199198043 —0.075757008 9.51E-08 
lichtenfelderi-catbaensis —0.338691314 —0.386903935 —0.290478693 3.37E-10 
luii-catbaensis —0.365692271 —0.425393393 —0.305991149 3.37E-10 
lichtenfelderi-huuliensis —0.201213788 —0.263738975 —0.138688602 3.37E-10 
luii-huuliensis —0.228214745 —0.299972435 —0.156457055 3.37E-10 
eye to ear distance (EE) 

huuliensis-catbaensis 0.08509 1637 0.059731157 0.110452117 3.37E-10 
lichtenfelderi-huuliensis —0.06878479 —0.094475902 —0.043093678 4.25E-10 
luii-huuliensis —0.101777804 —0.131262479 —0.072293128 3.37E-10 
luii-lichtenfelderi —0.032993014 —0.057865404 —0.008120623 0.00377334 
eye diameter (ED) 

huuliensis-catbaensis 0.085091637 0.059731157 0.110452117 3.37E-10 
lichtenfelderi-huuliensis —0.06878479 —0.094475902 —0.043093678 4.25E-10 
luii-huuliensis —0.101777804 —0.131262479 —0.072293128 3.37E-10 
luii-lichtenfelderi —0.032993014 —0.057865404 —0.008120623 0.00377334 
forelimb length (FLL) 

huuliensis-catbaensis —0.137477526 —0.199198043 —0.075757008 9.51E-08 
lichtenfelderi-catbaensis —0.338691314 —0.386903935 —0.290478693 3.37E-10 
luii-catbaensis —0.365692271 —0.425393393 —0.305991149 3.37E-10 
lichtenfelderi-huuliensis —0.201213788 —0.263738975 —0.138688602 3.37E-10 
luii-huuliensis —0.228214745 —0.299972435 —0.156457055 3.37E-10 
hindlimb length (HLL) 

huuliensis-catbaensis —0.137477526 —0.199198043 —0.075757008 9.51E-08 


The mito-nuclear data set of Liang et al. (2018) differed 
from all the above analyses in that their ML and BI anal- 
yses (79/0.99) placed the yingdeensis group as the sister 
group to a sister lineage comprised of the /uii group and 
lichtenfelderi group (87/1.00). 


Ancestral state reconstruction 


The AICc scores for the three transition rate models of 
the SCM analysis were ARD = 34.547134 and SYM and 
ER = 32.099451. The SCM analysis using either the SYM 
or ER model suggests that a karst habitat preference is 
the most probable ancestral condition for Goniurosaurus 
(57.0% probability), the kuroiwae group (62.7%), the 
luii group (90.0%), and the yingdeensis group (95.7%; 
Fig. 4). The probable ancestral condition for the /ichten- 
felderi group is non-karst (55.4%). The karst habitat pref- 
erence of G. kwanghua and G. zhoui of the lichtenfelderi 
group is considered to have evolved independently given 
that the ancestral condition of the /ichtenfelderi group 
and that of the most recent common ancestor of the sis- 
ter species G. lichtenfelderi and G. hainanensis was not 
karst-adapted (Fig. 4). 


Ecomorphology 


In both the species and habitat preference PCA analy- 
ses, PC1 accounted for 49.1% of the variation in the data 


set and loaded most heavily for limb length (FLL and 
HLL), snout length (SL), eye diameter (ED), interorbital 
distance (IO), head width (HW), and head length (HL). 
PC2 accounted for an additional 13.3% of the variation 
and loaded most heavily for body width (BW) and body 
height (BH) (Figs. 5, 6; Table 4). 

The PCA analysis of the karst-adapted Goniurosaurus 
catbaensis, G. huuliensis, and G. luii of the luii group 
demonstrates that their body shapes greatly overlap in 
morphospace despite there being several slight, but statis- 
tically significant mean differences among them (Fig. 5; 
Table 3). Additionally, none of the plots of the karst-adapt- 
ed species overlap with that of the granite stream-adapted 
species G. lichtenfelderi along the ordination of PC1. 

The PCA analysis using habitat preference as the de- 
pendent variable among the four species, showed that 
the karst-adapted and granite-stream-adapted species 
plot separately as before along taxonomic lines and that 
collectively, the former have significantly longer axil- 
la-groin lengths (AG); longer, wider, and thicker heads 
(HL, HW, and CH); longer snouts (SL); longer limbs 
(FLL and HLL); wider interorbital distances (IO); larger 
eyes (ED) and larger ear openings (AD) (Fig. 6). Many 
of these characters—longer head and snout, larger eyes, 
longer trunk, longer limbs—occur in many other dis- 
tantly related karst-adapted species of Cyrtodactylus 
(Grismer et al. 2016a, 2020b; Kaatz et al. 2021; Nielsen 
and Oliver 2017), indicating that these are convergent 
adaptions to a karstic life style within and between the 
gekkotan families. 
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Table 4. Summary statistics and principal component analysis scores for the mensural characters for Goniurosaurus catbaensis, G. huuliensis, G. luii, and G. lichtenfelderi. Abbreviations are listed in the Materials and methods. 
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0.12433 
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—0.11090 
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Discussion 


Geckos in general are particularly well-adapted to 
karstic landscapes (see Luu et al. 2016; Grismer et 
al. 2014, 2020a, 2021 and references therein; Google 
Scholar search using key words “karst” and “Gek- 
konidae”) and Goniurosaurus is no exception, being 
that 19 of its 23 species (83%) occupy karstic hab- 
itats (Grismer et al. 1994, 1999; Orlov et al. 2008; 
Ziegler et al. 2008; Yang and Chan 2015; Honda 
and Ota 2017; Zhou et al. 2018; Ngo et al. 2019a; 
Qi et al. 2020a, 2020b; Zhu et al. 2020a, 2020b). 
It is clear that karstic landscapes have played a sig- 
nificant role in the evolution and diversification of 
Goniurosaurus being that it is the probable ances- 
tral habitat preference for the genus and three of the 
four species groups. Even the ancestor of the non- 
karst adapted ancestor of the Jichtenfelderi group 
was karst-adapted (Fig. 4). Furthermore, within the 
species groups, the limited data herein would sug- 
gest that the karst-adapted species are specialized, 
range-restricted endemics (Grismer et al. 1994, 
1999; Orlov et al. 2008; Ziegler et al. 2008; Yang 
and Chan 2015; Honda and Ota 2017; Zhou et al. 
2018; Ngo et al. 2019a; Qi et al. 2020a, 2020b; Zhu 
et al. 2020a, 2020b). With the exception of G. licht- 
enfelderi, all the non-karst-adapted species are re- 
stricted to islands in the Ryukyu Archipelago (kuroi- 
wae group) or Hainan Island (ichtenfelderi group). 
It may be that the absence of competition and/or 
predators in these insular habitats widened the fun- 
damental niches of their ancestors and allowed some 
species to become more generalized in their habitat 
preference, which should be tested using new tech- 
niques combining phylogenetic history, character 
evolution, and ecological reconstruction programs. 


Systematics of the /uii group 


The ML and BI analyses of Liang et al. (2018) and 
the BEAST analysis herein (Fig. 4) recovered Go- 
niurosaurus catbaensis as the strongly supported 
sister species to the remainder of the /uii group. 
Whereas the ML and BI analysis herein, recov- 
ered G. catbaensis as the very poorly supported 
(60/0.51) sister species of the G. araneus plus G. 
gezhi clade (Fig. 3). Given that three of the five 
analyses strongly supported the former relationship 
and two analyses poorly supported the latter, we 
prefer the placement of G. catbaensis as the sister 
species to the remainder of the /uii group (Fig. 4). 
Given the very low nodal support of the latter, it 
essentially renders that portion of the tree a poly- 
tomy and as such, does not effectively contradict 
the strongly supported sister species position of G. 
catbaensis in the other trees. 

Goniurosaurus kadoorieorum of the luii group 
(represented by only 16S) is nested within G. /uii in 
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both the ML and BI analyses, rendering G. /uii polyphy- 
letic (Fig. 3). The same relationship was recovered in the 
16S phylogeny of Zhu et al. (2020a). This, and the lack of 
diagnostic characters separating G. kadoorieorum from 


G. luii (Yang and Chan 2015; Ngo et al. 2016), indicates 
the two species should be considered conspecific and as 
such, G. kadoorieorum Yang and Chan, 2015 is relegated 
here to a junior synonym of G. /uii Grismer, Viets, and 
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Boyle, 1999. In all analyses, G. huuliensis is consistently 
recovered as the sister species to G. Juii sensu lato and its 
species status is not questioned (Figs. 3, 4). 


Conservation 


Wide-ranging more inclusive studies pertaining to ecosys- 
tems management are becoming commonplace in light of 
climate change and widespread habitat destruction. Such 
studies reconcile data from a broad range of disciplines 
in order to address issues that may bear on ecosystems 
management. Foundational to many of these studies is a 
basic understanding of species ecology and habitat pref- 
erence—correlated here with ecomorphology (Cabral et 
al. 2009; Harfoot et al. 2014). Baseline information on 
habitat and microhabitat requirements of any species are 


paramount to understanding how they interact with, and 
navigate through, their environment (e.g. Grant and Grant 
2008; Greene 2005; Losos 2010) and as such, the contex- 
tualization of ecosystem management may ultimately turn 
on these simple points (Meiri 2018; Sinervo et al. 2010). 
Integrating the phylogenetic patterns of biodiversity 
and the morphological adaptations of habitat preference 
that, in part, underpin species radiations, can fundamen- 
tally contribute to conservation management programs 
(Grismer et al. 2020a, 2021; Erwin 1991; Vane-Wright 
et al. 1991; Williams et al. 1991; Vazquez and Gittleman 
1998; Moritz et al. 2000; Forest et al. 2007; Sgro et al. 
2010; Harvey et al. 2011; Rolland et al. 2012; Winter et 
al. 2012; Shaffer et al. 2015; Beaumont and Wang 2019; 
Fay et al. 2019; Holderegger et al. 2019)—especially in 
the karstic regions of northern Vietnam where anthropo- 
genic impact is degrading the habitat and reducing the 
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density of localized populations of Goniurosaurus (Ngo 
et al. 2019b). Northern Vietnam and many of its offshore 
islands in the Gulf of Tonkin, harbor large areas of frag- 
mented karstic habitats scattered across their landscapes 
(Cerrano et al. 2006; Do 2001, 2014; Luo et al. 2016; 
Ngo et al. 2019a) that are inhabited by an exceptionally 
large number of endemic plants and animals (Do 2001; 
Sterling et al. 2006; Clements et al. 2008; Luo et al. 2016; 
von Oheimb et al. 2017). The obligate restriction of many 
species to fragmented karstic environments—such as all 
species of the /uii and yingdeensis groups—functionally 
transforms these environments into habitat islands (Cle- 
ments et al. 2006, 2008; von Oheimb et al. 2017), which 
in some cases, bear an unprecedented degree of range-re- 
stricted endemism (e.g. Sgro et al 2012; Harvey et al. 
2011; Grismer et al. 2018a, 2021). 

Unfortunately, Goniurosaurus species are particularly 
attractive (Fig. 2) and over-harvested for the illegal pet 
trade (Stuart et al. 2006; Yang and Chan 2015; Ngo et 
al. 2019b). This is an additional threat to these range-re- 
stricted endemics from imperiled karstic environments 
(Grismer et al. 1997; Orlov et al. 2008; Ziegler et al. 
2008; Nakamura et al. 2014; Yang and Chan 2015; Hon- 
da and Ota 2017; Zhou et al. 2018; Ngo et al. 2019a; Qi et 
al. 2020a,b; Zhu et al. 2020a,b). In fact, in areas of China 
and Vietnam, many populations have suffered huge de- 
clines in numbers, or even extirpation at some localities, 
due to the illegal commercial pet trade (Stuart et al. 2006; 
Yang and Chan 2015; Ngo et al. 2019b, 2021). We hope 
that this study will bring more clarity to the plight of this 
genus and continue to serve ongoing conservation and 
management programs. 
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